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§ ■ Abstract 

We discuss various properties of the variational class of continuous matrix product states, a 
class of ansatz states for one-dimensional quantum fields that was recently introduced as the direct 

>. 

I/*-) _ continuum limit of the highly successful class of matrix product states. We discuss both attributes 

CO ' 

0^ ■ of the physical states, e.g. by showing in detail how to compute expectation values, as well as 

co ; 

properties intrinsic to the representation itself, such as the gauge freedom. We consider general 
£NJ \ translation non-invariant systems made of several particle species and derive certain regularity 

properties that need to be satisfied by the variational parameters. We also devote a section to the 
^ ■ translation invariant setting in the thermodynamic limit and show how continuous matrix product 

h: 

states possess an intrinsic ultraviolet cutoff. Finally, we introduce a new set of states which are 
tangent to the original set of continuous matrix product states. For the case of matrix product 
states, this construction has recently proven relevant in the development of new algorithms for 
studying time evolution and elementary excitations of quantum spin chains. We thus lay the 
foundation for similar developments for one-dimensional quantum fields. 
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I. INTRODUCTION 



Many revolutions and breakthroughs in quantum physics, and quantum many body 
physics in particular, were stimulated by guessing a suitable variational ansatz that cap- 
tures the relevant correlations for the systems under consideration. Feynman's ansatz 



for the roton in superffuid Helium 111. Il2| . the Bardeen-Cooper-Schrieffer wave function 
for superconductivity [| ] and the Laughlin wave function for the fractional quantum Hall 
effect 19|] are only a few prominent examples. For gapped one-dimensional quantum spin 



systems, the set of matrix product states 



10 



38 



is a very general ansatz that can 



describe a range of different phenomena and different physical phases, including normal 
symmetric and symmetry broken phases as well as the more exotic symmetry-protected 
topologically ordered phases such as the Haldane phase 



17 



18 



33] . Indeed, with the benefit 



of hindsight, we now understand White's powerful density matrix renormalization g roup 
algorithm |39L l40| as a variational optimization over the set of matrix product states 29|, |34|. 

Until recently, few equally general ansatzes that surpass mean field theory were available 
for extended quantum systems in the continuum, i. e. quantum fields. Numerical approaches 
require a finite number of degrees of freedom in order to fit the problem in the memory of a 
computer. For compact systems such as nuclei, atoms and molecules, an expansion in terms 
of a finite-dimensional basis is possible, but for extended systems this eventually results in 
a discretization to an effective lattice system. A new variational ansatz field theories in 

n 

d = 1 spatial dimensions was developed by Verstraete and Cirac in 2010 [37| . This ansatz is 
formulated in the continuum and does not require an underlying lattice approximation. It 
can be considered to be the continuum limit of a special subclass of matrix product states 
(MPS) and is therefore called the continuous matrix product state (cMPS) class. 

The aim of the current paper is to discuss in greater detail the properties of cMPS. 
Section UTI reviews the different definitions and representations of these states in the current 
literature. We then derive a set of regularity conditions that become relevant in the case of 
systems with multiple particle species in Section HTT1 Section |V] discusses how to (efficiently) 
evaluate expectation values with respect to these states. Section IVT1 is devoted to the gauge 
invariance and the existence of canonical forms in the continuous matrix product state 
representation for generic systems without translation invariance. We also discuss uniform 
continuous matrix product states in the thermodynamic limit and illustrate how continuous 



matrix product states possess a natural ultraviolet cutoff in Section |VII[ Finally, Section lVlIII 
provides an intuitive construction of tangent vectors to the variational set and discusses their 
representation properties as well, both for finite systems and in the thermodynamic limit. 
These tangent states are relevant when stud ying time evolution or elementary excitations 



along the lines of analogous MPS algorithms 



14 



16 



21 



3l| . We do not strive for absolute 



mathematical rigor, but merely attempt to explain in full detail the prerequisites for using 
cMPS in numerical algorithms. For example, due to the intrinsic difficulty of the various 
infinite-dimensional function spaces involved, we do not include a rigorous proof that the 
set of continuous matrix product states constitutes a smooth (complex) manifold and that 
the construction of a tangent space is justified. 



II. VARIOUS DEFINITIONS OF THE VARIATIONAL CLASS 
A. Setting 

Consider a quantum system defined on a one-dimensional continuum 1Z = [-L/2, +L/2] 
with length \1Z\ = L that accommodates q bosonic and/or fermionic particle species, which 
are labeled by the greek index a = l,...,q. Throughout this paper, we restrict to non- 
relativistic systems. A state of the quantum system containing N a particles of type a is 
then described by a square integrable function on n^i^te^; where t] a = +1 (—1) if 
particle species a is bosonic (fermionic) and IZ^^ (IZ^^) corresponds to the symmetric 
(antisymmetric) subspace of 7Z N , the Cartesian product of N copies of 7Z. The space of the 
square integrable functions on this domain is a Hilbert space that is denoted as 

Following the principles of second quantization, we now define the Fock space 

l4 F) = © ■ ■ ■ © H^-^ 1 -- (2) 

AT 1= o Nq=0 

which captures an arbitrary state of the quantum system. In addition, we denote the unique 
vacuum state as \Q) G e^ Q " 0}a=1 " " 9 . Particles of type a are created and annihilated at 
position x e 1Z with the operators i>' a (x) and i) a {x) with a = 1, . . .,q. These satisfy the 
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general commutation or anticommutation relations 

^a{x)^p{y) - ria^p{y)^ a {x) = 0, $ a (x)$p(y) - Va,^l(y)^ a (x) = S ajj3 6(x - y), (3) 

where rj a ^ = — 1 if both a and represent fermionic particles and r] a ^ = 1 when at least 
one of the two particles species a or j3 is bosonic. Clearly r\ a a = r\ a . We always write sums 
over the species index a explicitly and do not use Einstein's summation convention with 
respect to this index. 



B. Original definition 



A cMPS is defined to be the state 37| 



|^[Q, = tr ffiTexp 



+L/2 ^ 1 

dxQ{x) g) 1 + y2R a (x) ®i>t(x) 

Q=1 



|n>, (4) 



where CPexp is the path ordered exponential (that orders its argument from left to right for 
increasing values of x) and \Q) is the empty vacuum that is annihilated by ip a (x), Va = 
1, . . . , N. The trace operation acts on an auxiliary space C D , also called the ancilla space, 
where D is the bond dimension. The variational parameters correspond to the functions 
Q, R a : 1Z C DxD that take value in L(C Z) ) = C DxD , the space of linear operators acting 
on the ancilla space. For now, we do not impose any continuity or regularity conditions on 
these functions, and we refer to Section II I II for a detailed discussion. Finally, the boundary 
operator B G L(C D ) encodes the boundary conditions. For a system with periodic boundary 
conditions the boundary operator has full rank and is typically chosen to be B = Ip. In 
case of open boundary conditions, we can choose B = v^v^ with and «r D-dimensional 
boundary vectors. Note that the matrix functions Q and R a themselves need to satisfy 
certain boundary conditions which are imposed by the physical setting. We discuss this in 
more detail in Section IIVI 

More formally, we can identify the cMPS construction as a map between the function 
spaces TZ — > C DxD and the Fock space H^: 



:(TZ^C DxD ) q+1 



(F) . 

(5) 

(Q,Ri,...,R q ) i y |*[Q, 



The range of the map $ defines a variational set V c mps(d) C , where we often omit 
the explicit specification of the bond dimension. Henceforth, we compactly denote a cMPS 



\ty[Q, Ri, ■ ■ ■ , Rq\) as \^[Q, {R a }]). It will always be clear from the context how many and 
which particle species are present. The variational set V c mps(d) is not a vector space, since 
the representation of the sum of two elements \^[Q, {R a }}) + l^fQ', requires in the 

most general case a cMPS \^f[Q, {R a }]) G -^cMps(d) with bond dimension D = 2D, where 
we choose (Vx G [-L/2,+L/2]) 



Q(x) =Q(x)®Q'(x), 

R a (x) = R a (x) ®R' a (x), Va = l,...,g 

B = B © B'. 



The variational set does however contain almost complete rays of states, since for any state 
{Ra}}) £ V c mps(£>) an d any A G C = C \ {0} we can also represent A \^[Q, {R a }}) 
as a cMPS with bond dimension D as ^[Q', {R' a }]), where Q\x) = Q(x) + /i(x)lo and 



R' a (x) = R a (x) with 




exp / dx n(x) = A. 



A special case is obtained for A = 0, since this requires us to redefine Q(x) as Q'{x) = 
Q(x) — ocId. Hence, the null state is not contained within V c mps(d) but only in its closure. 
Correspondingly, the variational set V c mps(d') with D' < D is not a subset of V c mps(l>)- For 
example, if the boundary matrices are fixed to B' = l^i and B = (periodic boundary 
conditions), then a representation of the cMPS ^'[Q', {i?^}]) with bond dimension D' as 
a cMPS \$[Q, with bond dimension D > D' requires Q = Q' © (—00 x 1d~d') and 

R a = R' a © (0 x 11 £>_£)/), hence V c mps(d') is only included in the closure of V c mps(d)- Note 
that this differs from the case of MPS on the lattice, where Vmps(d') C Vmps(d) f° r D > D' . 
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C. Fock space embedding 



The embedding of \^[Q, {R a }]} £ V c mps(d) in the Fock space M.^ for finite \1Z\ can be 
made explicit by expanding the path ordered exponential as 

+00 „ 

MQ,{R a }]) = J2 dx!---dx N 

N=0 J -L/2<x 1 < - <x N <L/2 



tr 



B ( Q(xi) ® 1 + R c*i 00 ® ^ (xi 
«i=i 



qjv=1 



|n>. 



We can then expand the round brackets and reorder the sum in terms of the actual number 
of created particles by grouping subsequent occurrences of the Q term, so as to obtain 



+00 q „ 

ro,{iU]> = £ E / 

7V=0 ai,...,ajv=l 



dxi ■ ■ ■ dx 



N 



tr 



L/2<xi<-<x N <L/2 

BM Q (-L/2, x 1 )R 0ll (x 1 )M Q (x 1 ,x 2 ) ■ ■ ■ R aN (x N )M Q (x N , L/2) 

^M)^U x 2)---^ aN (x N )\n), (6) 



with 



+00 „ 

M Q (x,y) = J2 / dz 1 ---dz k Q(z 1 )---Q(z k ) = ?el'" Q ^ dz . 

Eq. shows how a cMPS can be interpreted as an superposition over the different particle 
number sectors in the Fock space. Note that this is not completely equivalent to the different 



sectors BL 



{N a } a= 

n 



l,...,q 



in the direct product construction of [Eq. (j2J)], since now only the 
total number of particles N = YH=i N a is fixed. If we define the iV-particle wave functions 
as 

( Xl ,...,x N ) = (n\i> ak (x k ) ■ ■ - ^(x^mQ^Ra}}) . (7) 
then we can infer from Eq. ([6]) that 



y oi,...,«jv 



a; at) 



tr 



BM Q {-L/2, x 1 )R ai {x 1 )M Q {x 1 ,x 2 ) ■ ■ ■ R aN (x N )M Q (x N , L/2) 



only when x\ < x 2 < • • • < xn- It can be extended to any other order of the arguments 
by reordering the annihilation operators in Eq. (J7J) according to the given commutation or 



anticommutation relations in Eq. ([3]). The non-relativistic kinetic energy requires that these 
functions are sufficiently regular, which together with the extension to arbitrary order of the 
arguments imposes certain non-trivial constraints on the matrix functions Q and R a that 
are to be discussed in Section IIHI 



D. The continuum limit of matrix product states 



The cMPS \^f[Q, {R a }}} was originally constructed in Ref. |37| as the continuum limit of 
a certain subset of MPS, where the subset was selected in such a way as to obtain a valid 
continuum limit. We explore this construction in greater detail and elaborate on some of 
the non-trivial implications regarding ultraviolet cutoffs and correlation lengths (infrared 
cutoffs) . 

We approximate the continuum 1Z = [-L/2, L/2] by a lattice £ with lattice spacing a 
and N = L/a sites, where we send a — > 0. On every site of the lattice we can create and 
annihilate particles of type a by acting with the creation and annihilation operators $ a {n) 
and c a {n). We can relate them to the field operators by 

/>(n+l)a 

c a {n) = / i a (x) dx (9) 

J na 

and its hermitian conjugate. The local basis on site n thus consists of the states |0) (no 
particles), \a) n = c£ (n) |0) n , \a,(3) n = c£ (n)ct(n) |0) n , . . . On this lattice, we can define an 
MPS with matrices A s (n) where s can take values 0, a, (a, /3), ... If the local basis is 

infinite-dimensional, this MPS definition is only formal, i.e. it cannot be used for practical 
computations. In the limit a — > 0, the number of sites L/a in the lattice £ goes to infinity. 

On an infinite number of lattice sites, two arbitrary MPS are generally orthogonal due 
to the (infrared) orthogonality catastrophe [3j. Since we now aim to create quantum field 
states within the Fock space HL^ , we need to restrict to a special subset of MPS where the 
total number of particles is finite (on average, so that (N) is finite). Since a finite number 
of particles has to be distributed over a diverging number of sites L/a, most of the sites in 
the lattice £ are empty on average. So A has to be the dominant matrix, and it turns out 
that the cMPS \V[Q,{R Q }]) E can be obtained from the continuum limit (a — > 0) of 



the MPS \^[A}) G He by identifying ^ a {na) = & a {n)/y/a and 

A°(n) = t D + aQ(na), 
A a {n) = yfaR a (na), 

A^\n) = { 2 (10) 
^R a (na) 2 , a = (3 



together with \Q) = |0) = <S>„ e £ |0) n , Vn = —L/2a, —L/2a + 1, . . . , +L/2a — 1. This equiv- 
alence can be obtained from a Taylor expansion of the exp-operator, although this is only 
completely rigorous when the entries of Q and R a are finite and the operators i]A(x) are 
bounded {i.e. not for bosons). Most results for cMPS in the remainder of this chapter can be 
derived from this correspondence with MPS, but we attempt to derive these results directly 
in the continuum as much as possible. 

The correspondence with MPS is useful for concluding that the entanglement of one half of 
the chain with the other half (in the case of open boundary conditions) is limited by the upper 
bound log-D. By restricting to MPS within a single Fock space in the thermodynamic limit, 
we avoid the orthogonality catastrophe. The infrared orthogonality catastrophe of MPS in 
the thermodynamic limit would turn into an ultraviolet catastrophe when this infinitely- 
sized lattice L would correspond to the continuum limit of a finitely sized continuum 1Z. 
Physically, the ultraviolet catastrophe is avoided because the finite number of particles 
induce a physical cutoff a p h ys that is given, not by the lattice spacing a — > but by a p h ys = p^ 1 



with p = (N) / L the particle density |24l|. The presence of a physical length scale can be 
detected from the physical dimensions of Q and R a , which are given by [Q] = i~ x and 
[R] = £ -1 / 2 with £ a generic length dimension. The nature of the physical cutoff a p h ys 
and its relation to Q and R a is discussed in Section IVIII for the translation invariant case, 
where it can unambiguously be defined. Shifting the cutoff from the lattice spacing a to 
a physical value a p h ys is a very important step in the definition of cMPS. MPS with finite 
bond dimension D have a finite amount of entanglement to which corresponds in general a 
finite range of fluctuations £/a, where £ denotes the correlation length. Hence, they have 
in general a finite dimensionless correlation length £ = As a is scaled to zero while £ 
remains finite, the physical correlation length £ would also scale to zero. It is because the 



physical cutoff is shifted to a finite value a p h ys (with thus a p h ys /a. — > oo) that cMPS are able 
to combine a finite amount of entanglement with a finite physical correlation length £ (with 
thus £/a — > oo but with £/a p h ys finite). The physical correlation length £ is also computed 
in Section IVHI for the translation invariant case. 



E. Alternative construction through continuous measurement 

Rather than trying to construct a cMPS as the continuum limit of a MPS, we could also 
try to directly define the continuum limit of the processes that define MPS. Unfortunately, 
the process of sequential Schmidt decompositions has no straightforward generalization to 
the continuum and neither has the definition of valence bond solids. One can however define 



a continuum version of the sequential generation process that creates MPS 35|, based on 
the paradigm of continuous measurement {(J. The resulting process for creating cMPS is 



described in Ref. 



28 



I and is here summarised for the sake of completeness. 

As in the discrete case, let the ancilla start in a state G EI anc iii a = C D . This ancilla 
can be interpreted as a resonating cavity with D internal levels, in which there is a particle 
source that creates particles of type a (a = 1, . . . , q). These particles gradually leave the 
cavity due to cavity losses. Since particles leaving the cavity at different times occupy 
different positions in space at a given time (since they travel at a certain speed which we 
set equal to one), the resulting configuration of particles can be interpreted as a static 
spatially distributed quantum state. For a compact cavity (i.e. a zero-dimensional system), 
the resulting quantum state is one-dimensional. As an abstraction of this physical process, 
a (d — l)-dimensional cavity can be used to encode a d- dimensional holographic quantum 



28 



for the general case, and henceforth restrict to the d — 1 case that 



state. We refer to Ref. 
produces cMPS. 

Between two particle emissions, the cavity evolves according to a Hamiltonian K G L(C D ) 
(a Hermitean D x D matrix), whereas the physical state outside the cavity does not evolve. 
By observing the particles that are emitted from the cavity, we are continuously measuring 
the state of the cavity (i.e. ancilla). The state of the cavity at time t is encoded in the 
particle distribution at position x = —t. It was shown that the resulting configuration of 

10 



particles outside the cavity is given by 

u^exp -i / dxK(x)®i + y2iR a (x)0i)l{x)-iR a (x) i ®ij; a (x))v K \n), (11) 

V J ~ L /* a=l J 

where the ancilla is projected onto the state ul at the end of the measurement, in order to 
decouple it from the physical state. The resulting expression does not yet correspond exactly 
to Eq. (jl]) but it can easily be brought in the required form by using the Baker- Campbell- 
Hausdorff formula on every infinitesimal patch of the path ordered exponential. We then 
obtain that the state in Eq. flTTJ is contained within V c mps ; as it is equal to \^[Q, {R a }]) 
for the specific choice 

1 N 

Q(x) = -iK(x) --J2 R <*(rfRa{x). (12) 

1 a=l 

We recall that K(x) is a Hermitian matrix. Generic cMPS can be brought into this form by 
using the gauge invariance of the cMPS representation, as discussed in Section [VTl 
This construction allows us to introduce a unitary operator U (y, z) G L(C D <g> M) 

U(y,z) = ?exp i^' l j z ^K(x)®l + ^\R a (x)®^l(x) -m a (x) ] ®i) a {x^ . (13) 

Being a unitary operator, it conserves the norm of i>r <S> \Q). This does not imply that the 
cMPS \^[Q, {Ra}}) with Q given by Eq. ( Fl2|) is automatically normalized to unity, because 
the definition also involves a projection to v^. But the unitarity of U(y, z) in Eq. ffl3|) does 
guarantee that \^[Q, {Ra}]) can easily be normalized and has no norm that diverges or goes 
to zero in the large volume limit. 

From a physical perspective, this construction is important as it clearly sketches the holo- 
graphic properties of the cMPS. The physical state of a one- dimensional system is described 
by a zero-dimensional boundary theory. The spatial coordinate of the physical system acts 
as a time coordinate in the boundary theory. The physical state is created because the 
boundary theory interacts with the physical system, where the position of the interaction 
shifts linearly in time. This interaction results in the boundary theory not being at equilib- 
rium. Instead, the boundary theory is subject to dissipative dynamics, as will become clear 
in the following section. This holographic property is of course strongly related with the 
intrinsic area law for entanglement that is present in cMPS. 
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F. Path integral representation 



Recently, it has also been illustrated that we can break up the path ordered exponential 
in the definition of |*[Q, {««}]> and insert rations of the identity in order to obtain a 
path integral description of the same state[5|. The easiest way to insert an identity is by 
first introducing a second quantized version of the ancilla by making the substitution 

Q(x) i y Q(x) = Q h \x)b% R a (x) ^ R a (x) = W a k {x)b% (14) 

with bj and b^ annihilation and creation operators for bosonic or fermionic particles in level 
j — 1, . . . , D of the ancilla. The resolution of the identity can now be expressed in terms of 
coherent states. However, the ancilla Hilbert space is now an infinite-dimensional Fock space, 
whereas the original ancilla space was only C D and corresponds to the single-particle sector 
of this Fock space. Because the operators Q(x) and R a {x) are particle- number preserving 
with respect to the ancilla, we can restrict the whole path integral to the single particle sector 
by either choosing appropriate boundary conditions. If \u) denotes the ancilla zero-particle 
state, then a restriction to the single particle state is obtained by identifying 

B = B j ' k b] \u) (u\ b k . (15) 

If we introduce the coherent states 

10) =exp M ( 16 ) 

then we can write the identity as 

1 r D 

J 3=1 

Following the standard recipe, we can then obtain the path integral description of \^[Q, {R a }}) 
as 

MQ,{R a }]) = 

D0£>0* (0(+L/2) f J B0(-L/2)) e" 



x exp 



L/2 



Q 

J2(<t>Kx)R«(x)<t>(x))fia(x)}dx |n>, (18) 

a=l 
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where <fi{x) is a D- dimensional vector function with components <f)j{x), j = 1, . . . , D. This 
path integral representation can serve as a useful starting point for generalizations of the 
cMPS, e.g. by replacing the second quantized auxiliary system by a true field theory, so that 



23 



If this field theory 



this becomes the cMPS analogon of the construction in Ref. Q and 
is a conformal field theory, it is then very close in spirit to some model states for Quantum 
Hall Systems Q]. 



III. REGULARITY CONDITIONS 

In Eq. ((7j) we have defined the iV-particle wave functions 0a ll ... J a JV (xi, . . . , Xn). For x\ < 
■■■ < x^ these are completely specified by Eq. (jSj). However, for general choices of the 
matrix functions Q and R a , the extension of Eq. (jSj) to all orders of its arguments does 
not automatically satisfy the required properties that a physical iV-particle wave function 
should satisfy. For example, the iV-particle wave functions should be different iable in each 
of its arguments if the state has to produce a finite non-relativistic kinetic energy. 

However, there is no need to work with the Fock space expansion of Eq. §6$). We can 
check the regularity of the iV-particle wave functions by immediately evaluating the kinetic 
energy in second quantization. For further reference, we first define 



U(x, y) = CPexp 



ry ( 1 
J &z lQ(z) ®t + ^2R a (z) ®ft a {z) 



(19) 



where U(x,y) 6 L(H®C D ) with C D the ancilla space, i.e. it is a D x D matrix of operators. 
Unlike the operator U(y, z) defined in Subsection III CI the operator in Eq. ffl9|) is not unitary. 
It only equals the unitary version when acting on and if Q(z) is given by Eq. (Tl2|) . In 
addition, we define a closely related set of operators U a (x,y) (a — 1, ... ,q) as 



U a (x,y) = ?exp 



dz < Q(z) ® 1 



13=1 



(20) 



In order to compute any expectation value, which is the topic of the next section, we need 
to be able to act with the field annihilation operators ip a {x) on the state \^[Q, {R a }}). If 
we are able to drag tp a ( x ) through the path-ordered exponential, it then acts on which 
is annihilated by any field operator. We can now use Eq. (lAlip as derived in Appendix lAl 
where B = ^(x), Ai(z) contains both Q(z) ® 1 and any term Rp{z) (g> ^(z) for which 
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Va,p = 1, and A 2 (z) contains the terms Rp(z) ® ^(z) for which r} a) p = — 1. We then obtain 

$ a {x)U(-L/2,+L/2) - U a (-L/2,+L/2)^ a (x) = U a {-L/2,x)R a U{x,+L/2) 
which immediately results in 



ipa(x) \ V[Q,{Rp}]) = tr BU a (-L/2, x)R a (x)U(x, +L/2) \Q) . 



(21) 



Hence, acting with an annihilation operator of type a at position x not only lowers a matrix 
R a {x), but also transforms the path ordered exponential U(—L/2,x) into U a (—L/2,x), 
because we had to take the particle statistics into account for bringing ij) a {x) to the position 
where it could lower R a (x). 

The non-relativistic kinetic energy operator T is given by 

-+L/2 



t(x) dx, 



L/2 



(22) 



where the kinetic energy density t(x) at position x is given by 



N 1 



a=l 



2m n \ dx 



[x 



dx 



[x 



(23) 



Hence, a finite kinetic energy expectation value ($[Q, {R a }]\T\^[Q, {R a }]) requires that the 
state ^-(x) \V[Q, {R a }}) has a finite norm. Differentiating Eq. (I2ip and using Eq. ( IA2I) . we 
obtain 



—Mx)MQ,{Rp}}) 



tr 



dR a 



BV a (-L/2,x)[ [Q{x),R a {x)] + -^-(x) )U(x,+L/2) 



|fl> 



tr 



B7 a (-L/2,i)(^ [^^(x)i? 



/9=1 



®^J(x) )U(x,+L/2) 



|n> . (24) 



The term on the first line can be shown to have a finite norm (see next section), provided 
of course that R a {x) is a different iable function with a well-behaved derivative dR a (x)/dx 
at any x & 1Z. Since the term on the second line of Eq. (124)) has particles of any species 
(3 = l,...,q being created at the fixed position x, this term is not normalizable. Put 
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differently, || (d^(x) / dx) \^[Q, {R a }}) || 2 contains a divergent contribution 5(0) (in position 
space), unless we impose the regularity condition 

r) a> pRp{x)R a {x) - R a {x)Rp(x) = 0, V16 11. (25) 

Hence the matrices R a should have the same statistics as the particle creation oper- 
ators to which they couple. For systems with a single species of bosons, the condi- 
tion in Eq. (1251) is automatically fulfilled. For systems with multiple species of bosons, 
it requires that any two matrices R a (x) and Rp(x) at the same spatial point x com- 
mute. If a is a fermionic particle species, the corresponding matrix R a {x) has to satisfy 
Ra{x) 2 = 0, Va; G 1Z. When two particles of fermionic type a approach each other, there 
is a corresponding factor R a (y)'J } exp(J^ dx Q(x))R a (z) in the iV-particle wave function 
(x±, . . . , y, z, . . . , x N ). For y -» z, the exponential factor continuously evolves 
towards Id, so that the fc-particle wave function continuously becomes zero. Hence, the 
finiteness of the kinetic energy requires that two fermionic particles of the same type cannot 
come arbitrarily close together and thus imposes Pauli's principle. 

Differentiability of the wave function is sufficient for a finite kinetic energy, which is by 
far the most important physical requirement of the wave function. We can also impose 
higher regularity constraints on the iV-particle wave functions. Since these do in general 
not arise from physical considerations, we postpone this discussion to Appendix [B] While 
the resulting conditions are interesting from an algebraic point of view, they are in general 
hard to satisfy with finite-dimensional matrices. For practical applications, satisfying the 
original condition in Eq. (125]) . as imposed by the finiteness of the kinetic energy, should be 
sufficient. 

We conclude this subsection by investigating what else can be learned from the physical 
considerations concerning particle statistics. The regularity conditions [Eq. (125|) ] already 
require that the matrices R a behave as the corresponding operators ip a in terms of commu- 
tation and anticommutation relations. In a physical system, we should not have fermionic 
condensates, i.e. (^[^(x)^) = if particle species a is fermionic. This is a consequence 
of the invariance of an physical Hamiltonian H under the action of the parity operator P, 
which flips the sign of any fermionic operator (Pip a (x)P = r) a>a ip a (x)) and is thus idempotent 
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(P = P 1 = pt). We can construct P as 



exp 



a fcrmionic 



cxp 



E 



q fcrmionic 



7v 



(26) 



Physical states satisfy P \ ^>) = e 1 ^ \ ^>) , where the idempotence of P requires <fi = or = tt. 
Physical states thus consist completely of a superposition of states, all of which have either 
an even or an odd number of fermions. Imposing this same property for cMPS requires 
one to explicitly incorporate the Z2 symmetry (with group elements {I, P}) in the matrix 
structure of R a and Q. Since P \^/[Q, {R a }]) = \^[Q, {r]a,aR a }]), we should also be able 
to define a virtual operator P G L(C Z) ) such that PQP^ 1 = Q and PR a P~ l = r] ayCe R a . 
This operator can in principle be ^-dependent, but we should then be able to apply a local 
gauge transformation (see Section IVT)) in order to make P space- independent. In addition, 
it is clear from the definition that P is idempotent (P = P _1 ). If we can assume that 
P is diagonalizable, then P divides the ancilla space C D into a sector with positive parity 
(eigenspace of eigenvalue +1) and a sector with negative parity (eigenspace of —1). A global 
gauge transformation brings P into the diagonal form 



(27) 



1d(+) d(+ ) xD (-) 
Od(-)xd(+) -Id(-) 

with P/ + ) + = D. The required transformation behavior of Q and R a then imposes 
the following decomposition 



Q 



R« 



Rn 





Od(+)x£>(-) 


°D(-)xD(+) 




d(+) 


Od(+)xd(-) 


°D(-)xD(-) 




Od(+)x£>(+) 


d(+") 


d(-+) 


Od(-)xd(-) 



(particle species a is bosonic), 



(particle species a is fermionic). 



(28) 
(29) 
(30) 



In the cMPS \^f[Q, {R a }]), all contributions with either an even or an odd number of fermions 
in Eq. ([6]) drop out, depending on the boundary matrices B. If only states with an even 
number of fermions are allowed, B should have a decomposition as 

B {+) D (+) XjD (-) 
0d(-)xd(+) 



B 



(31) 
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whereas a decomposition of the form 



B 



0, 



(32) 



is required to select only states with an odd number of fermions. 

IV. BOUNDARY CONDITIONS 

We have already mentioned in Section [TT] that the type of boundary conditions — open 
or periodic — is encoded in the rank of the boundary matrix B. For a system with periodic 
boundary conditions, B has full rank and is typically chosen to be the identity (B = Id)- 
Since periodic boundary conditions identify the points x = —L/2 and x = +L/2, it is natural 
to assume that the matrix functions Q and R a are also single- valued, i.e. Q(—L/2) = 
Q(+L/2) and R a (-L/2) = R a (+L/2) for all a = 1, . . . , q. 

For a system with open boundary conditions, it is suitable to work with a boundary 
matrix of the form B = v^v^, i.e. the rank of B is one. However, in the case of open 
boundary conditions physical requirements impose additional conditions on the TV-particle 
wave functions of Eq. (JSJ). Typically, a finite system is interpreted as being embedded in an 
infinite system and having an infinitely strong potential energy outside of the interval TZ, 
i.e. v(x) = +oo for x < —L/2 and x > +L/2. The single particle wave functions that build 
up the Fock space are zero outside 1Z. A finite kinetic energy imposes continuity, and thus 
requires that the single particle wave functions are zero at x — ±L/2. Consequently, the 
resulting A-particle wave functions have to produce zero as soon as one of the arguments 
Xi is equal to ±L/2. Since this has to be true for any configuration of the remaining N — 1 
particles, we obtain that we have to impose 

v[R{-L/2) = and R{+L/2)v R = 0. (33) 

where a partial differ- 



A more detailed discussion of these conditions is presented in Ref. 
ential equation for the evolution of Q and R a under real or imaginary time dynamics is de- 
rived. In order to solve this partial differential equation, it needs to be complemented by the 
proper boundary conditions as given above. Throughout the remainder of this manuscript, 
we assume that we are working with cMPS where the matrix functions Q and R a satisfy 
the required conditions. 
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We now also have to discuss whether we can completely fix the boundary matrix B, 
or whether its entries should be included within the set of variational parameters. While 
B = Ip represents a fixed choice that is well-suited for the case of periodic boundary 
conditions, we will see in Section I VI I that it is beneficial to include one of both boundary 
vectors v l or in the set of variational parameters in the case of open boundary conditions. 
In order to have a uniform notation, we do not explicitly denote this dependence in the 
notation for the state \^?[Q, {R a }])- Note that it is impossible to absorb the boundary 
vectors into the matrices Q(—L/2), R a (—L/2) and Q(L/2), R a (L/2) in the case of open 
boundary conditions. More generally, unlike in the case of generic MPS on finite lattices, it 
is for cMPS impossible to use a space-dependent bond dimension D(x), since the required 
continuity of D in combination with its discrete character enforces a constant value. 



V. COMPUTATION OF EXPECTATION VALUES 



This section is concerned with the computation of expectation values of normally ordered 
operators. We have already illustrated how to act with annihilation operators and derivatives 
thereof in the Section HTT1 With a MPS, the computation of expectation values boils down to 
a contraction of the physical indices in the network. In the continuum, however, the intuitive 
notion of physical indices is a bit lost. We therefore start by computing the overlap of two 
cMPS \^[Q, {R a }}), \^[Q', {R'a}}), which are given as an expansion in Fock space [Eq. (JBJ)]. 
It is clear that the basis states ip^ (x\) ■ ■ ■ ipa N ( x N) \&) are automatically orthogonal for 
different N, and further that 



WavM ■ ■■ ^1(2/1)^(^1) ■ • • 4>l N ( x N)\ty = 

<Wi ■ ■ • 8 aN ,i3 N S(xi - yi) ■ ■ ■ 6(x N - y N ), (34) 
due to the ordering of the arguments X\ < ■ ■ ■ < xn and y\ < ■ ■ ■ < yN- We thus obtain 

+00 q „ 

mQ',{R' a }]MQ,{R a }]) = Yl E / da*.- da* 

N=0 {a u ...,a N }=l J-L/2<xi<-<x N <+L/2 



tr 



/ nxi \ / P+L/2 

BVexpl Q(z)dz i4 1 (zi)---i4 JV (zjv)5 , exp / Q(z) dz 

\J-L/2 J \Jx N 



x tr 



•I'l 



B7ex P [ I Q'(z)dz) R' ai (x 1 )---R' aN (x N )Texp 

L/2 



-L/2 



Q(z) dz 



x N 
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Using trivial direct product identities such as ti[A] tr[B] = tr[A <g) B], (AB) <g) (CD) = 
(A ®B)(C® D) and exp(A) ®exp( J B) = exp(A <g> t D + t D <g> B) for D x D matrices A, B, 
C and D, the previous expression can be rewritten as 



+00 H „ 

mQ',{R' a }]MQ,{R a }]) = J2 E / d Xl ---dx N 

N=0 {au ^ aif}=1 J-L/2<x l <...<x N <+L/2 



tr 



;5 <g> J B)Texp ( / [Q(z) ® 1 D + 1 D ® Q'(»] dz (^(xx) ® 

-i/2 

+L/2 



(Ra N (xN) ® i?' (x Ar ))0 5 exp 



[Q(z) g) 1 + 1 ® Q'(z)] dz 



•''A' 



Reverting the expansion of the path ordered exponential that lead to Eq. ([6]), results in 



'MQ'AKmiQAR*}}) 



tr 



(B®B)9explJ lQ(x)®t D + l D ®Q>(x) + Y,RM®R'a(x)]dx 



(35) 



-V2 a =l 

From the expression above, we can deduce that in the computation of expectation values 
(Q' = Q, R' a = R a ) a central role is played by the local transfer matrix T(x) defined as 

N 

T(x) = Q(x) ® 1 D + Id ® + ^ ^0*0 ® ^0*0- ( 36 ) 

Q!=l 

To this transfer matrix, we can also associate linear maps 7^ : L(C D ) t— > L(C D ) and 
: L(C D ) h-> L(C D ) that map virtual operators f (D x D matrices) to 

N 

0<*)(/) = Q(x)/ + /Q(z)+ + E^(^)/^(^) t » (37) 



N 



7 {x) (f) = fQ{x) + Q(a:)V + E R^fRaix). 



(38) 



The transfer matrix T(^) is of course strongly related to the transfer matrix E(n) = 
^2 s A s (n) ® y4 s (n) that features in expectation values with respect to MPS on the lattice. 
Indeed, if |^[^4]) is the MPS with matrices A as in Eq. (ITUl) . then the transfer operator T(x) 
is related to the transfer operator E(n) of the MPS \V[A]) by E(n) = 1 + aT(na) + 0(a 2 ). 

The expectation value of any normally ordered operator O =: : can now be 

computed by first acting with all annihilation operators ip a (x) on the ket \*&[Q, {Rp}}) as we 
did in the Section [TTT| and similarly acting with the creation operators on the bra. The result 
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of this is the insertion of some operators acting on the virtual system at the corresponding 
positions, with operators U(x,y), U a (x,y) or U^^x^y) connecting them. The expectation 
value is obtained by "contracting the physical indices" , which results in the inserted virtual 



25[, whereas the 



operators in the ket combining with those in the bra at the same position 
contraction of the part in between the local insertions result in a path ordered exponential 
of the transfer matrix. However, to incorporate the particle statistics, we also need to define 
generalized transfer operators as 

N 

T a (x) = Q(x) ®t D + t D ® Q{x) + ^2rj a , f} Rf)(x) ® R p {x) 



/3=1 
N 



7=1 



(39) 
(40) 



Note that T r 



x 



T(x) since if a a = 1. Given this recipe we can, for example, evaluate 



the correlation function 



G°^(x,y) = (nQ^R a }Ui(x)My)mQ,{R a }}) 



8(x — y) tr 



+ 9(y — x) tr 



(B ®Bpe f -^ T ^ {z)dz (Rp(y) <g> t D )9e f y Ta(z)dz 

x {t D ®Rjx)p^ LI2T ^ dz 

(B ® B) 7e J -^ TaAz) dz (t D ® R^(x)) T " w dz 

x (Rp(y)®l D pe£ L/aT W Az 



(41) 



All quantities in this expression, if we could store and manipulate variables with a fully 
continuous x-dependence, are D 2 x D 2 matrices. Since such matrices need to be multiplied, 
this is an operation with computational complexity of 0(-D 6 ), or 0(D 5 ) if we exploit the 
tensor-product structure. 

For physical systems, we can further simplify Eq. f l4Tj) . When only bosonic parti- 
cle species are present, all r\ a ^ = 1 and T = T a = T a ^. In case of the presence of 
fermionic particle species, we should incorporate the Z 2 parity symmetry discussed in the 
Section IIHI We can then define an idempotent parity superoperator P = P <g> P and we 
obtain PTP = T, as well as PT a P = T a and PT Qj/3 P = T a> p. This allows to conclude that 
(tf(Q, {R a }]\^l{x)^ p {y)\^[Q, {R a }]) = whenever the particle species a and f3 have differ- 
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ent statistics. When a and (3 are both bosonic or both fermionic, it is clear that T Qj( g = T 
and T a = Tg. 

In the case of open boundary conditions, we can define virtual density matrices l(x), r(x) G 
~h(C D ) which are defined through the initial conditions /(—L/2) = v^v^ and r(+L/2) = 
VrV^ and the first order differential equations 



^-l(x)=7^(l(x)), 



and 



-r(i) = -TW(r(i)). (42) 
ax v ' 



To these density matrices l(x) and r(x) we associate vectors \l(x)), \r(x)) G C D <8> C-° in the 
ancilla product space. Formally, the solution is given by 

(l(x)\ = (l(-L/2)\?e I -^ J{y)dy , 
\r(x)) = ?e^ +L/2ny)dy \r(+L/2)). 



We can then write 

mQ,{R a }]MQ,{R a }})= \l(-L/2) 



CPexp 



-L/2 



L/2 



T(x) dx 



r(+L/2) 



(l(x)\r(x)) = tr [l(x)r(x)] , Vx <E7Z. 



(43) 



From the correspondence with completely positive maps, it can be shown that the solu- 
tion l{x) and r{x) of Eq. (142 p starting from positive definite initial conditions /(—L/2) and 
r(+L/2) are positive for any i6K (see Theorem 3 in Ref. |2Cj). The norm is thus guaran- 
teed to be positive. Note that, for the special parameterization of Q(x) in the continuous 
measurement interpretation [Eq. (1121)]. we can write the determining differential equation 
for r(x) as 



dx 



r( x ) = -7 {x) (r(x)) 



^ N N 

i[K(x),r(x)] - - J2{MrfR*(x),r(x)} + ^ R a {x)r{x)R a {x)^ . (44) 



Q=l 



a=l 



This is a master equation in Lindblad form 20( describing the non-equilibrium Markov 
dynamics of the ancilla (i.e. the cavity). Starting from a pure state r(L/2) = v-rV^ at 
t = —x = —L/2, it evolves through interaction with the physical system (via the interaction 
operators R a ). At a general time t = —x, the density matrix r(x) is no longer pure: non- 
equilibrium evolution is a dissipative process. Note that the evolution is trace preserving, 
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since tracing the equation above results in dtr[r(x)]/dx = 0. In addition, the corresponding 
map 7& satisfies 7^(1 D ) = 0. 

In systems which only contain bosons, all r] a ^ = 1 and there is no need to introduce 
T a (x), T ai p(x), etc. As an alternative to the general recipe described above, we can then 
also deduce all expectation values of normally ordered operators O =: 0[{t/>£ }, : from 

a generating functional Z[{ J a }, { J a }} as (see Ref. |28|) 



o 



5 



8Jp) [5 J/3 

with 5 /SJ a the functional derivative with respect to J a , and 



Z[{Ja},{J«}} 



(45) 



Z[{J a },{J a }]=ti 



r+L/2 

(B®B)?exp\ J dxT(x) 

N 



+ Ja{x)[Ra{x) ® Id] + J a (x)[l D <g> R a (x)} 

a=l 

which for a system with open boundary conditions results in 
Z[{J a },{J a }]= (l{-L/2) 

+ ^2j a (x)[R a (x) ® 1 D ] + J a (x)[l D ® R a (x)} 



(46) 



+L/2 

?exp{ I dxT(x) 

L/2 



N 



Let us now illustrate this approach 
system with open boundary conditions 

H =f +V+W = 

-+L/2 x / d 



r(+L/2) . (47) 



dv defining a generic Hamiltonian for a single-boson 



L/2 



dx I — ibHx) ) I —ibix 

2m Vdx v ' Vdx 



-L/2 



L/2 



dxv(x)^ (x)tjj(x) 



^ r+L/2 r+L/2 

+ -/ dx dyw(x,y)4)\x)4) ] (y)%i)(y)%i)(x) (48) 

1 J -L/2 J -L/2 

describing particles with mass m that interact with an external potential v(x) and with each 
other through two-particle interaction w(x,y). 

Using Eq. (1451) we find (henceforth omitting the arguments Q and R in the state |\l/)) 



($\^(x)ip{x)\y) = (l(x)\R(x) ® R{x)\r(x)), 



(49) 
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and 



9{y - x) (l(x)\R(x) ® R^jVe^ dzT ^R{y) ® R^y)\r{y)) 

+ 0(x - y)(l(y)\R(y) ®Rjy)?e J y dzT{z) R(x) ®R(x)\r(x)). (50) 

Denning R x l \x) = R(x)U(x)R(x) for every x G [— L/2, +L/2] and solving 

^(J#(v)l = (J4°(v)TO (51) 

for every y G [x, L/2], we can write the expectation value of the potential and interaction 
energy as 

+L/2 



(y\V\V) = dxv(x)(l(x)\R(x)®R(x)\r(x)), (52) 

J -L/2 

= /^dx r L \yw{x,y){Rf{y)\R{y)®R^y)\r{y)). (53) 



J-L/2 Jx 

To evaluate the expectation value of the kinetic energy, we compute 



\dx J \dx J x-*v axay 

= lim 



>y dxdy 



lim — 

da; 



9(y - x)(l(x)\(l D ® i?(x))ye^^ T W(^) ® lij)|r(2/)) 

+ 6»(a: - y)(l(y)\(R(y) ® l D )!Pe^ T(2) (l D ® ^))}|r(z)) 

0(y - x) (Z(x) | (1 D ® 15)) OW* d * T(2) 

[¥(</), ® 1 D ] + {dR(y)/dy ® 1 D ) )|r(y)) 
+ 0(x - y) (Z(y) 1 1 [T, ® 1 D ] + (dR{y)/dy ® lj,) 

x ?e^ d2T(z) (l D ®^)|r(x)) 
We have used the defining equations [Eq. (142]) ] in the computation of d(l(y)\/dy = (l(y)\T(y) 



x 



and d\r(y))/dy = -T(y)\r(y)). Since T(y) = Q(y) ® 1 D + 1 D ® Q(y) + ® we 
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obtain [T(y),R(y) ® 1 D ] = ® Id and thus 



lim 



9(y - x)(l(x)\l D <8> ([Q(z),i2(x)] + di?(x)/dx)lPe^ d2T(2) 

x ([Q(y), i?(y)] + dR(y)/dy) ® l^r^/)) 
+ 0(a! - </)(%) |([Q(y), fl(y)] + di2(y)/dy) ® l D Te^ d2T « 



x1d®(1 d ® [Q(x), R(x)] + dR(x)/dx) \r(x)) 

where we used the same trick. Note that derivatives with respect to the Heaviside functions 
(which would produce a diverging contribution 5(x — y)) nicely cancel for both derivatives 
with respect to y and to x. As noted in the Section IIII[ the regularity condition Eq. ( 12 5 p is 
automatically fulfilled for the case of a single boson. We thus obtain 

-L/2 



mm = — / dx{l{x)\{[Q(x),R(x)}+dR(x)/dx) 
2m J_ L/2 



® ([Q(x),R(x)] + dR(x)/dx)\r(x)). (54) 

Note that this result could also be obtained by the general strategy outlined at the beginning 
of this section, i.e. by acting directly on the cMPS with the operators ip(x) and d^(a;)/dx 
and only afterwards computing the expectation values. However, the generating function 
approach is very general and relates nicely to the standard approach that is used to compute 
expectation values in quantum field theory. As for the definition of the state itself, we can 
also write the generating functional using a path integral, which can be useful for analytic 
computations or Monte Carlo based evaluation strategies. 



VI. GAUGE INVARIANCE 



As with a MPS, the map ^ associating a physical state \^[Q, {R a }\) G to the matrix 
functions Q : TZ — > C DxD and R a : 1Z — > <C DxD is not injective, i.e. the representation is 



led gauge invariance was rigorously discussed in terms 



not unique. For MPS, this so-ca 
of principal fibre bundles in Ref. [15j. Such a rigorous treatment for the case of cMPS is 
severely complicated by the fact that both the domain and the codomain of the map ^ are 
now infinite dimensional. Therefore, it is beyond the scope of the current manuscript, as 
noted in the introduction. We thus proceed in an intuitive way. 
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We do expect the existence of a local gauge transformation g : 1Z — > GL(D, C), i.e. a 
position-dependent invertible matrix g(x), that acts on the matrices Q(x) and R a {x) while 
leaving the physical state \^[Q, {R a }) invariant. While it is hard to extract the correct trans- 
formation formulas for Q and R a from the original cMPS definition in Eq. (j3J), people with a 
background in Yang-Mills gauge theories might recognise Q as the connection that generates 
parallel transport by comparing the iV-particle wave functions of the Fock space embedding 
[Eq. ([S])] to Wilson lines with insertions of charges transforming according to the adjoint 
representation, or from recognizing the action of the path integral formulation [Eq. ( [18]) ] as 
a Yang-Mills action with a covariant derivative ^- + Q(x). The gauge transformation for a 
cMPS is thus given by 



Q(x) = g(x) 1 Q(x)g(x) + g(x) 



dx 



x), 



R(x)=g~\x)R(x)g(x) 



(55) 



While we prefer the continuum derivation, these transformation formulas can also be ob- 
tained by using the correspondence with MPS [Eq. (fTUl) ] and the well-known gauge trans- 
formations for MPS 15] 

A°(n) = g{{n - l)a)- 1 A°(n)g(na) 

= g((n — l)a) _1 g(ria) + ag({n — l)a)~ 1 Q(na)g(na) 
dg~ l 



t D + a 



na)g{na) + g{na) 1 Q{na)g{na) 



A a (n) 



dx 

g((n-l)ay 1 A a (n)g(na) 
V^g(na) _1 ii! Q (na)g(na) + 0(a 3 ^ 2 ) 
) = g{{n-l)a)- 1 A^g{na) 



+ 0(a 2 ), 



| [R a (najRpina) + rj a ^Rp{na) R a {na)} + 0(a 2 
1R a {na) 2 + 0(a 2 ), 



/3 

a = (3 



Indeed, using dg~ 1 (x)/dxg(x) = —g~ 1 (x)dg(x)/dx, we reproduce the transformation for- 
mulas of Eq. f[5"5"j) . To have an invariant physical state \^[Q, {R a }}) = \^[Q, {Ra}}), we 
also need to transform the boundary matrix as B = g(L/2)~ 1 Bg(—L/2). When B is 
fixed, we need to restrict to gauge transformations that satisfy the boundary condition 
g(L/2)~ 1 Bg(-L/2) = B (e.g. g(L/2) = g(—L/2) for B = 1 D ). In addition, we also re- 
quire the function g : 1Z — >• GL(D, C) to be second order differentiable in order to have 
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new matrix functions Q(x) and R a (x) which have a well-defined first order derivative. The 
regularity condition of Eq. (1231) is not modified by the gauge transformation and puts no 
further constraints on the set of allowed gauge transformations. Since this condition follows 
from physical considerations which are left invariant by gauge transformations, it would be 
strange if we obtained a different result. 

As for MPS, we can use the gauge fixing conditions to impose a certain canonical form on 
the matrices Q{x) and R a (x). Suppose we want to impose a gauge fixing condition such that 
Q(x) is of the form in Eq. ({12"]) . corresponding to the cMPS construction from continuous 
measurement. It is equivalent to the left orthonormalization condition of MPS and boils 
down to imposing 

i 

Q{x) + Q(x) t + R a {x)^R a {x) = 

a=l 

for every x e TZ. Inserting the explicit form of Q(x) and R a (x) in terms of the original 
Q(x), R a (x) and g(x) [Eq. (1531) ]. we obtain that g{x) should be a solution of the differential 
equation 

= (g-\x)y g ~ 1 (x)Q(x)+Q(xy (g-\x)) ] g-\x) 

+ Y,R a (x^ {g~\x)) ] g~\x)R a (x) (56) 

a=l 

= 7^ [(g-\x)y g^x) . 

Clearly, this differential equation only determines g(x) up to a unitary prefactor. Put dif- 
ferently, for any solution g(x) of this equation, g'(x) = u(x)g(x) with u(x) a unitary matrix 
is an equally valid solution. We can use the remaining gauge freedom u{x) G U(D) to 
diagonalize r(x) at every point x, hence obtaining the left- canonical form. 

However, at this point it becomes important to discuss the boundary conditions that 
should be satisfied by solutions g(x). If the boundary matrix B is fixed, we need to impose 
g~ 1 (+L/2)Bg(—L/2) = B. This is a highly non-trivial condition and it is not certain that 
such solutions exist. For periodic boundary conditions with B = 1^, it logically results in 
g(+L/2) = g(—L/2). Translation-invariant states with periodic boundary conditions can be 
subjected to the the same treatment as the translation-invariant states in the thermodynamic 
limit, which are discussed in the next section. Henceforth, we restrict to the case of open 
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{g-\x)) ] g'\x 



boundary conditions with B = VrV^. From this, we can derive the conditions 

vlg{-L/2) = av{ g-\+L/2)v R = -v R 

a 

for some non-zero a G C. However, we can easily fix a = 1 by substituting g(x) <c— g'{x) = 
g(x)/a, since the constant gauge transformation alp acts trivially on Q and R, i.e. it is 
within the kernel of the gauge group action. Nevertheless, the resulting boundary condi- 
tions are still highly non-trivial and it is not assured by the standard theory of differential 
equations that there exist solutions satisfying both conditions simultaneously. Hence, it is 
better to restrict to a single boundary condition such as g(—L/2) = to and do not im- 
pose any condition on g(+L/2). The value of g(+L/2) is then completely determined by 
the differential equation (up to the unitary prefactor). Consequently, we then also have to 
transform the right boundary vector as v R = g~ 1 (+L/2)v R . This implies that v R is part 
of the variational degrees of freedom, and should also be included in e.g. the variational 
optimization for finding ground states. Note that the boundary conditions for g(x) are in- 
herently imposed by the representation of the state, and are not related to or influenced by 
the physical conditions that need to be satisfied by Q and R, as discussed in Section IIV1 

Alternatively, we can also impose the right orthonormalization condition, which boils 
down to 

N 

Q{x) +Q{x) ] + ^R a {x)R a {x) ] = (57) 



a=l 



and implies that 



Q(x) = -iK(x) -\Yj R a (x)R a (x)^ (58) 

a=l 

with K(x) a Hermitian matrix. Starting from an arbitrary cMPS with matrices Q(x) and 
R a (x), we obtain new matrices Q(x) and R a ( x ) according to Eq. fl55l) . which satisfy the 
above relations if g(x) is a solution of 

d q 

— [g(x)g(x) j ] = -Q(x)g(x)g(x) ] - g(x)g(x) ] Q(x) ] - R a (x)g(x)g(x) ] R a (x) ] 

dx (59) 

= [g(x)g(xy] . 

Clearly, for any solution g(x), we obtain a family of solutions g'(x) = g(x)u(x) with u(x) G 
U(-D). This unitary freedom can be fixed by diagonalizing l(x), resulting in the right- 
canonical form. As for the left-canonical form, one has to pay careful attention to the 

27 



boundary conditions that need to be satisfied by g. For a system with open boundary 
conditions, the easiest solution is again to include one of the boundary vectors in the set of 
the variational parameters and also transform it under the action of the gauge transform. 

Note that we can also define a gauge transformation g(x) for the cMPS \^f[Q, {R a }]) £ 
A^ c mps so that 



Q(x) = g(x) 1 Q(x)g(x) + g{x) 



It is sufficient to choose 



g{x) = CPexp 



+L/2 



Q(y)dy 



go 



(60) 



(61) 



with g some arbitrary integration factor that is fixed by the boundary conditions. For 

and we also 



f^Q(y)dy 



example, if we require g(—L/2) — Id then g = (^CPexp 
need to transform <— = g(+L/2)~ 1 vn = g^v-R. Hence, the cMPS can now be written 

as 



M{R a }}) = <CPexp 



r+L/2 N 



v R \n) 



(62) 



This formulation is close in spirit to the bosonic mean field ansatz 

-+L/2 



\(p) = exp 



-L/2 



tp{x)ip\x) dx ] \Vt) 



with ip a scalar (complex-valued) function, since it identifies the mean field ansatz with a 
cMPS with bond dimension D — 1. This mean field ansatz lies at the basis of the Gross- 
Pitaevskii equation 



13 



32|, that is still used today with great success. All variational 



degrees of freedom are now contained in the matrices R a {x) (and -Dr), and all gauge degrees 
of freedom have been eliminated. However, we do not employ this particular choice of 
gauge in the remainder of this manuscript as it also has some downsides. For example, 
translation- invariant states \^f[Q, R a ]) can be obtained by choosing the matrices Q and R a 
x- independent (see next subsection). However, this particular gauge transformation maps 
the x-independent matrices R a to x-dependent matrices R a (x) = e + ® x R a e~ Qx , so that 
translation invariance is less easily recognized. 

VII. TRANSLATION INVARIANCE AND THE THERMODYNAMIC LIMIT 

When using cMPS to approximate ground states of translation invariant Hamiltonians, 
we can restrict to the subclass of uniform cMPS \^(Q, {R a })), which are obtained from 
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taking Q(x) = Q and R Q (x) = R a constant x-independent D x D matrices in \^[Q, {R a }}). 
This approach is valid both for a finite system with periodic boundary conditions (B = 
or for a system in the thermodynamic limit (\7Z\ = L — > oo or thus 1Z — > R), where 
the precise value of the boundary matrix B should be irrelevant and should not appear in 
any normalised expectation value. We henceforth restrict to the latter case. The trans- 
fer operator T = Q ® Id + Id ® Q + Yla=i ^ a ® ^ a a ^ so becomes translation invariant 
and 7 exp [ f dx T] = exp[T(2; — y)\. The normalization of the state \ty(Q,R)) is given by 
lim^ootr [(B ® B) exp(TL)] . If /x = max Aeo .( T ){3ft(A)}, where cr(T) denotes the spectrum 
of T and 3? the real part, then (&(Q, {R a })\^(Q, {R a })) ~ lini£_ s . 00 exp(/xL). Normalizing 
this state by multiplying it with exp(— \xV) results in Q ^— Q — fi/21 D and T ^— T — /il, so 
that the new transfer operator T has at least one eigenvalue for which the real part is zero 
and no eigenvalue has a positive real part. Let us assume that the eigenvalue A with 9ftA = 
is unique. If \r) is the corresponding right eigenvector, then we can write the eigenvalue 
equation as T(r) = Ar with r the associated virtual density matrix. Hermitian conjugation 
learns that 7(r>) = Xr\ so that the uniqueness of the eigenvalue with 3ftA = implies that 
A = A = and = e 1? V, where we can choose the phase of the eigenvector so that r is 
Hermitian. Similarly, the virtual density matrix / associated to the left eigenvector \l) can 
also be chosen Hermitian. 

Having a unique eigenvalue zero and 3?(A) < for all other eigenvalues A corresponds 
to the generic case, as can be better appreciated by referring to the well-known results 



for MPS 



10 



15 



obtained by identifying[2 



30(. Indeed, a full categorisation of the eigenvalue structure of T can be 



T = lim-lnE (63) 

a->-0 a 



with E the corresponding transfer operator of the uniform MPS ^(A)) with A related to Q 
and R a as in Eq. ffTUj) . The set of MPS with a well-defined thermodynamic limit correspond 
to the injective or pure MPS, for which the transfer operator E has a single eigenvalue 1 
that maps to the eigenvalue zero of T. The corresponding left and right eigenvectors (7| and 
|r) correspond to strictly positive Hermitian operators I and r (i.e. they have full rank). All 
other eigenvalues of E lie strictly within the unit circle and map to eigenvalues of T with 
strictly negative real part. If the left and right eigenvectors corresponding to eigenvalue 
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are normalized such that (l\r) = 1, then lirn^oo exp(TL) 



r)(l\ and we obtain 



(*(Q,{R a })\*(Q,{R a })) = (l\B®B\r). 



(64) 



In expectation values of local operators, this overall factor always appears, but the rest of the 
expression will not depend on B. Hence, the independence is cancelled by considering nor- 
malized expectation values, or by henceforth choosing B such that ($f(Q, {R a })\^f(Q, {R a })) = 
(l\B®B\r) = 1. 

For uniform cMPS, the gauge invariance is restricted to global transformations Q <— Q = 
qQq 1 an d R a <— Rot — gRag' 1 with g e GL(C, D). This gauge transformation can be 
used to impose the left or right orthonormalization conditions. Left orthonormalization 
boils down to fixing the left eigenvector I of eigenvalue to / = 1^, which results in 
Q = —\K — 1/2 X]«=i B) a R a with K a Hermitian matrix. The remaining unitary gauge 
freedom can be used to diagonalize r, bringing Q and R a in the left-canonical form. The 
right-canonical form is obtained analogously. In principle, an exact computation of the left 
and right eigenvectors / and r corresponding to the eigenvalue with largest real part A of 
the transfer operator T are computationally costly operations [0(-D 6 )]. By using an explicit 
parameterization of the left-canonical form in terms of R a and the Hermitian matrix K, we 
know exactly that A = and I — tjj. It is then possible to obtain r with an iterative solver 
with computational efficiency 0(D 3 ). 

By imposing the physical requirements discussed at the end of Section IHIl we can define 
the parity superoperator P as in Section [V] Since FTP = T, we can expect that the left 
and right eigenvectors \l) and \r) corresponding to the zero eigenvalue satisfy (Z|P = (/| and 
P|r) = \r), or thus PUP = I and PrP^ = r. Note that we can always choose the gauge such 
that P is Hermitian. In addition, it is easy to prove that T a also has an eigenvalue zero 
even if a refers to a fermionic particle species so that T a ^ T. The corresponding left and 
right eigenvectors are in that case given by l a = IP = PH and r a = Pr = rP\ whereas they 
equal / and r if a is a bosonic particle. 

We can now evaluate correlation functions as 



c a , p {x,y) = (y(Q,{R a }Mi(x)My)\*(QAR*})) 

= 6{x - y)(l\[Rp ® l D ]e ra(x - v) [lD®TQ\r) 

+ 6(y - x)(l\[l D ® 7Qe r -to- x \R p (g> l D ]\r), (65) 
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where we have used the physical requirement T Q>y g = T and T a = Tp for non-vanishing corre- 
lation functions (see Section |V|). The correlation function C a ^(x,y) is translation invariant 
and we define C a ^(x, y) = C a ^(y — x). When a is bosonic and f3 fermionic, we automatically 
have C a ^(x) = if the parity considerations from Section HTT1 are correctly built in. In the 
long-range limit, we obtain lim^i^oo C a> p(x) = (l\Rp <8> tD\ r a)(la\^D R a \r). When both a 
and /3 refer to fermionic particle species, this limiting value is automatically zero (also under 
the assumption that parity is correctly built into the matrices). When both indices refer to 
bosonic particles, a non-zero value is possible in the case of Bose-Einstein condensation. We 
should then define a connected correlation function C a ^(x), which decays exponentially as 
limui^oo C a fi{x) = 0(exp[— |x|/£ c ]) with £ c = (SftAi)" 1 , where Ai is the eigenvalue of T a with 
second largest real part (i.e. skipping eigenvalue Ao = 0). Clearly, C a> p(x) is continuous 
at x = 0. We can then compute the first derivative, which is only continuous at x — if 
we impose the regularity conditions in Eq. (|25|) . This is another way to derive these condi- 
tions. If Eq. ( 1251) is satisfied, then the second derivative of C a ^(x) at x = (which gives 
the expectation value of the kinetic energy density t up to a factor — l/2m) is finite and 
automatically continuous. The third derivative is then finite but will not be continuous in 
general, without imposing further conditions as discussed in Appendix [Bl 
We define the Fourier transformed correlation function 



In order to evaluate n a ,p(p), it is important to separate exp(T a x) into two parts. The first 
part is given by E> a = \r a )(l a \, the projector onto the eigenspace corresponding to eigenvalue 
of T Q , and yields a singular contribution to the integral. If we define the complementary 
projector Q a — 1 — S a , then the remaining part 



is well behaved in the Fourier transform, since all of its eigenvalues decay exponentially x. 
If we then introduce the notation Q a (— T a ± ip)~ 1 Q a = (—T a ± ip) p , which is well defined 




(66) 



with 




(67) 



exp(T Q x) -S a = Q Q exp(T Q x)Q Q = Q a exp(Q a T a Q Q a;)Q ( 
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even at p = because the zero eigensector of T a is projected out, we can rewrite n a< p(p) as 



n a ,p(p) = 2ir5(p)(l\t D ® R a \r a )(l a \Rp ® t D \r) 

+ {l\[t D ® i4](-T Q + rp) p [i^ ® l D ]|r) 

+ (/|[^ ® 1 D ](-T a - ip) p [l D ® i£]|r). (68) 

The first term is only present for bosonic particles that have condensed. It would also 
disappear in the Fourier transformation of the connected correlation function C(x, y). If we 
define Fourier transformed field operators &(p) — no confusion between the state |^) and 
the momentum-space operator \P should arise — as 

i r+oo 

$(p) = -= / dx^(x)e- ipx , (69) 

then it is easy to see why we have used the suggestive notation n a $ for the Fourier transform 
of C a> p. We obtain 

^(Q,{R a })\^l(p)^{p f )mQ,{R a })) = 5{p-p')n a , p {p). (70) 

Hence, n a ^{p) describes the occupation number of momentum levels. The large-p behavior 
of n a fi{p) follows from the regularity of C a ,p{x). At first sight, Eq. (l68~l) might seem to 
decay as Q(p~ l ). However, if the regularity conditions in Eq. ( 1251) are satisfied, then the 
momentum occupation number n a ^{p) has to decay as 0(p~ A ) for large values of p. We can 
show this explicitly. For \p\ larger than the eigenvalue of T a with the largest absolute value, 
we can expand (— T a ± ip) p as 

(-T. ± *)' = Ti— £ = =F^ + ^ ± 4 - ^ + O(P-). (71) 

P V P J P P P P 

We now have to show that by plugging this expansion into Eq. (168|) . the first three terms 
vanish. The first term is trivial, if particle type a is bosonic so that Q Q = t — \r)(l\. For the 
fermionic case, one has to employ the parity conservation. Using the regularity conditions 
of Eq. (|25p and t] a ^ = r]p n for non-vanishing correlation functions — a and /3 are of both 
bosonic or both fermionic — we can show that 

T a [R p ® l D ]\r) = [Rf,®lD]T\r) + [Q,Rf,]®l D \r) = [Q,Rp] <8> t D \r) 
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and similarly 



T a [l D ®TQ\r) = l D ®fQ,TQ\r), 
(l\[Rp®l D ]T a = (l\[Rp,Q}®t D , 
(l\[l D ®TQT a = (l\t D ®[R^,Q\. 



These results can be used to show that both the second and third term in the expansion 
vanish when they are plugged into Eq. (|68|) . The first no n- vanishing term is thus of order 
p~ 4 . Because n a ,p(p) is a dimensionless quantity, this asymptotic behavior allows us to 
introduce a momentum cutoff A as 



A 4 = lim|Aa,/?(p)| = \(l\[t D ®R a ]Tl[R^t D ]\r) + (l\[R^l D ]Tl[t D ^R a ]\r)l (72) 



where the absolute value is not required if we use j3 = a. The eigenvalue spectrum of T a 
thus provides a definition for an ultraviolet cutoff scale a = A -1 . Rather than defining the 
ultraviolet cutoff scale a = A -1 through the total particle density 



we have now defined a UV cutoff scale A based on the large momentum behavior of the 
momentum occupation number n a ^{p). 

For two pure uniform cMPS \^{Q, {R a })) and \^(Q', {R' a })) we can define a superop- 



erator T mixed = Q' <g> t D + 1 D <g> Q + J2%=i K ® R<* so that the (V(Q, {R a })\V(Q', {R' a })) 



decays as limL^ +00 exp(AL), with A the eigenvalue with largest real part of T m ; xe d- If the 
two uniform cMPS are inequivalent, 9R(A) < and there is an infrared orthogonality catas- 
trophe. If 3?(A) = 0, then we can define a phase <fi = 3(A) and a gauge transformation 
g G GL(Z);C) such that Q' = gQg^ 1 + i(f> and R' a = gR a g~ x . With / being the right 
eigenvector corresponding to eigenvalue A = i<p of T mixe d, g can be obtained as g = fr' 1 . 

Let us also illustrate how to compute the expectation value of a translation invariant 
Hamiltonian. The generic Hamiltonian in Eq. (I48p becomes translation invariant for v(x) = 
v and w(x,y) = w(y — x) with w(x) = w(—x). Since the uniform cMPS is extensive, 
expectation values are proportional to the volume and it makes more sense to compute the 
expectation values of the kinetic, potential and interaction energy densities t, v and w. We 



p— ¥00 




(73) 
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obtain 

(V(Q,{R a })\i\V(Q,{R a })) = ^(!\[Q,R] ® \QMr), (74) 
(y(Q,{R a })\v\*(Q,{R a })) =v(l\R®R\r), (75) 

r+oo 

(V(Q,{R a })\w\V(Q,{R a })) = / ^^(^(/li?®^ 2 ^®^^). (76) 

If 10(2;) has a Laplace transform £[u>](<7) = / + °° dzw(z) exp(— az) that is defined for Sftcr > 0, 
we obtain 

(y\w\V) = (Z|i2®S£H(-T)i2®S|r). (77) 

Note that translation invariance has allowed the parameterization of a field with a continuous 
number of degrees of freedom by a discrete number of degrees of freedom. Having / and r , the 
computational cost is 0(D 6 ) when long-range interactions are present, since we then have to 
compute an arbitrary function L [w] of the transfer operator T, unless w is such that there is 
an exact or approximate (iterative) strategy for evaluating the action of £ [w] (— T) on a vector 
efficiently. One particular example is the case of strictly local interactions w(x—y) ~ 5(x—y). 
The interaction energy (density) can then be computed with computational complexity of 
0(D 3 ) just like the potential and the kinetic energy density. 



VIII. TANGENT VECTORS OF CONTINUOUS MATRIX PRODUCT STATES 



A. Generic case 



For MPS, a new algorithm for time evolution and variational optimization (via imaginary 
time evolution) was recently constructed using the time-dependent variational principlefl^. 
An essential ingredient of this algorithm is the study of (infinitesimally) small variations of 
MPS, i.e. the set of MPS tangent vectors. Indeed, it was rigorously proven that the set 
of MPS can be given the structure of a variational manifold with a well-defined tangent 



space [15j by eliminating some singular points or regions. While we do expect the same theo- 
rems to hold for cMPS, the infinite dimensionality of the parameter space and Hilbert space 
might require a different proof strategy, especially in the absence of translation invariance. 
As noted several times before, this would be beyond the scope of this paper. 
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Given the practical use of tangent vectors, we nevertheless proceed, albeit in a more 
intuitive manner. Let us assume that we do have an open subset of cMPS with fixed 
bond dimension D that constitute a (complex) manifold -M c mps C H. At any base point 
{R a }]) £ -Mcmps; we can construct a (holomorphic) tangent space 7|^[Q { flce }])A^ C MPs C 
H. If the collective index i = 1, . . . , D 2 is used to combine both virtual (matrix) indices 
(a, (3) and we use the summation convention with respect to this index, a general tangent 
vector |$[V, {W a }\ Q, {R a }]) in T mQt{Ra}]) M cM ps can be defined as 

mv,{Wa};Q,{R*}}) = \^ [Q ' {Ra}] [v,{w a }}) 



l/2 \ v '5<?{x) tt 6R mJ (78) 

|n>. 



+L/2 

' BU(-L/2, x) I V{x) Qi + ^Wpix)® tpl(x) ) U(x, L/2) 



dx tr 

L/2 



Because of the gauge invariance discussed in Section IVIt not all variations in Q and 
R a result in changes of the physical state. Consequently, not all linearly independent 
choices of the matrix functions V and W a result in linearly independent tangent vectors 
\$[V, {W a }; Q, {R a }]). Let Q(rj) and R a (rj) (Va = 1, . . . , q) be a one-parameter family of ma- 
trix functions, so that Q(rj) : 1Z i— > C DxD : x h- > Q(x; rf) and similarly for R a {rf). If we define 
Q(0) = Q : x ^ Q(x), Ra(0) = R a : x ^ # a (x) together with dQ/dr/(0) = V : x i-> V(ar) 
and di? a /d?7(0) = W a : x H )■ W Q (:r), then we can write 



±-MQ(r,),R a ( V )]) 



= mV,{W a };Q,{R a }]). (79) 

T7=0 

If we now choose a one-parameter family of gauge equivalent states, so that Q(x\rj) = 
g(x;r])~ 1 Q(x)g(x;r]) + g(x, r/) -1 ^*'^ and R(x;r]) = g(x;r])~ 1 R(x)g(x;ri), where the one- 
parameter family of gauge transforms is given by g(x; rj) = exjp(rjh(x)) and h(x) G fjt(C, -D) = 
C DxZ) , Vx G 72., then we can use the gauge invariance of the cMPS representation to obtain 
\^[Q(x;rj),R(x;rj)]) = \V[Q(x),R(x)]) and thus 

where the maps M^' and DSfjf$ (Va = 1, . . . , N) are given by 

M [ S ] [h](x) = [Q( x ),h(x)] + ^(x), ^ ] [h](x) = [R a (x),h(x]]. (81) 
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The maps M.® and thus establish a linear homomorphism from functions h : 1Z — > 

qI(C,D) = C DxD to the kernel of the representation |$[V, {W a }; Q, {R a }]} of the tangent 
space \^[Q, {R a }]) € 2]*[Q,{j? a ,}])A / f C MPS- Put differently, the representation of cMPS tangent 
vectors has a gauge invariance under the additive transformation law V <— V + [h] 

rn 1 

and W a <— W a + N a | [/i] . In all of the above, we have considered B fixed. The gauge 
transformation g(x) then has to satisfy the boundary condition g(+L/2)Bg(—L/2)~ 1 = B, 
which also imposes a boundary condition on the set of allowed functions h(x), namely 

h(+L/2)B - Bh(-L/2) = 0. (82) 

In particular, for periodic boundary conditions with B — Ip, we obtain that the generator 
h : 1Z — > gl(D, C) should satisfy periodic boundary conditions h(+L/2) = h(—L/2). 

We now restrict to the case of open boundary conditions and discard the explicit reference 
to the base point \^f[Q, {Ra}]) in the notation of tangent vectors. To take full advantage 
of the gauge freedom, we noted in Section [VI] that is better to include one of the boundary 
vectors in the set of variational parameters. We thus generalize our definition of tangent 
vectors by also including variations with respect to e.g. the right boundary vector v&. We 
write 

Mv,{w a },w R }) 

= W R ■ V VR \V[Q, {R a }}) 

p+L/2 / r N r \ 



0=1 



-L/2 

= vlU(-L/2,+L/2)w R \Q) 

r+L/2 ( N \ 

+ / dx vlU(-L/2, x) V(x) ® 1 + V Wf)(x) ® ipUx) U(x, L/2)v R \Q) . 

J -L/2 V 0=1 J 

(83) 

Let us revisit the gauge freedom for the new tangent vectors of Eq. (j83p . The state 
\$[V, {W a }, %]) is invariant under the additive gauge transformation V ^— V + M$[/i], 
W a <— W a + Na,<s> [h] and w R ^— w R + m$ [h] with 

m$[h] = -h{+L/2)v K . (84) 

Since is still fixed, the gauge transformation has to satisfy the boundary condition 
g(—L/2) = Id, so that its generator h(x) satisfies h(—L/2) = 0. 
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The overlap between two tangent vectors is given by 

($[F, {W a },m]mv, = wii(L/2) w ' R 

-+L/2 1 

dx {I (x) I V W' a (x) <8> W a (x) \r (x) ) 
V2 a= i 

+L/2 /-+L/2 9 

dx / dy (l(x)\[V'(x) ® 1 D + ^CW ®^R]Te^ d2T ^ 



L/2 



+L/2 



a=l 



[1 D (g) + £ i^fo) ® W^)] |r(j/)) 



55) 



-L/2 



dx / dy{l{y)\[l D ®V{y) + Y J My)®WM\^ AzT 



L/2 



a=l 



[V (x) ® 1 D + ^ W^(x) ® ^j] |r(x)) . 



a=l 



It defines a metric for the manifold -M c mps and features in any coordinate-invariant ex- 
pression involving cMPS tangent vectors. We can use the gauge freedom in the repre- 
sentation of tangent vectors to simplify the expression above significantly. The count- 
ing argument for the gauge degrees of freedom is now less rigorous as in the discrete 
case. In general, we have D 2 parameters in h(x) to eliminate D 2 degrees of freedom from 
{V(x), W\(x), . . . , W q (x)} at every point x. However, this is only correct if all linearly in- 
dependent algebra-valued functions h : 1Z — > g((C, D) map to linearly independent matrix 
functions [M,f\ {NjJ* 1 }]. Let us show that by substituting V{x) <- V(x) = V(x)+Wi*[h](x) 
and W a (x) W a (x) = W a (x) + N Q) $[/i](x) (Va = 1, . . . , q), we can indeed impose D 2 con- 
ditions, such as the left gauge fixing condition: 

N 

(Z(x)| V{x) <g> t D + Waix) ® R a (x) =0. (86) 



n=l 



This requires that h is a solution of 
d 



dx 



[l{x)h(x)] =7 {x) [l(x)h(x)} - 



l(x)V(x) + j2R a (x)H(x)W a 



x) 



a=l 



(87) 



which together with the boundary condition h(—L/2) = results in the solution 



(l(x)h(x)\ 



L/2 



dy(l(y)\ 



v(y) ® i D + w <*(y) ® R ^y) 



a=l 







Texp 


/ T(z)dz 




Jy 
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This equation gives a solution for l(x)h(x). We can extract h(x) by multiplying with /(x)^ 1 
to the left. The left density matrix l(x) should be positive definite and hence invertible 
for every x > —L/2. However, at x — —L/2 it equals l(—L/2) = v^v^ and thus becomes 
singular. Nevertheless, the limit \im x ^_L/2 h(x) should be well defined since the right hand 
side of the equation above, which is being multiplied with /i(x) -1 , will have a similar scaling. 
Alternatively, we can also impose a right gauge fixing condition 



N 



V(x) ® 1 D + W a {x) ® R a (x) 



a=l 



\r(x)) = 0. 



(89) 



Finally, we remark that the tangent space T\^Q^ Rct ^J^i cMPS spanned by the states of 
Eq. 083]) contains the original cMPS \W[Q, {R a }]), e.g. by choosing V = 1/L, W a = and 
w R = or by choosing V = W a = and ujr = vr. Both choices are related by a gauge 
transform with h{x) = (x/L + 1/2)1 D . For a general tangent vector \&[V, {W a }, iwr]), we 
obtain 



'MQ,{Ra}}MV,{W a },w R }) = vll(L/2)w R 

-L/2 



+ / dx{l(x)\V{x)®l D + VW" Q 

J ~L/2 a=1 



[x) ® R a (x)\r(x)). 



(90) 



If we fix the gauge according to either the left or right gauge fixing prescription, the 
second term cancels. We can restrict to the orthogonal complement of \^[Q, {R a }]) in 
T\<i,[Q t {R a }])M C MPS, which is denoted as T^^^n^M^ps, by further imposing 



v R l(L/2)w R = 0. 



(91) 



B. Uniform case 

We specialize again to the case of translation invariant systems in the thermodynamic 
limit. While the parameter space is now finite dimensional, it is fruitful to still consider 
the full tangent space to the manifold of all (translation non-invariant) cMPS at the special 
uniform point \^(Q, {R a })). This boils down to allowing space-dependent matrix functions 
V(x) and W a (x) in the definition of the tangent vectors. We can then decompose the 
full tangent space into sectors T$ p of momentum p e K by introducing Fourier modes 
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V(x) = Ve ipx and W a (x) = W a e ipx , resulting in 



\%(V,{W a };Q,{R a })) = (V,{W a })) = 

/ dxe ipx vlU{-oo,x) lv®l + ^2w a ®i>l{x) \ U{x, +oo)v R \Q) . (92) 

Note that the boundary vectors "Ul,r are irrelevant for the bulk properties of these states, and 
they are therefore not included in the set of variational parameters in the thermodynamic 
limit. Consequently, we also do not need to differentiate with respect to one of them in 
order to define the tangent space. 

We can also compute the overlap between two of these tangent vectors and obtain 

9 



/+oo " 
&xjW- p) *{l\Y,W' a ®W a \r) 

/+oo r+oo 1 
dx / dye i{p ' x - py) (l\[V ®l D + J2 W L®FQe {v ~ x)T 
-OO J X i 



a=l 

x [l D ®V + ^2R a ®\fiQ\r) 

a=l 

/+oo PX 1 

dx / dye i( - p ' y - px \l\[l D ®V + J2 R <*® VQe {x ~ y)T 
-oo J -co a=1 

Q 

[V® l D + ^W' a ®TQ\r). 



x 

a=l 



If we again resort to the decomposition of Eq. fIVIIj) . we can further evaluate this to 
2n5(p'-p) Ul\^2w^W a \r) 



a=l 



+ (/| [V ® 1 D + W a ® R a ](-T + fp) p [l D ®V + J> a ®W^r) m 

a=l a=l 
1 1 

+ (l\[l D ®V + Y^R a ®W^{-T-i P ) p [V' ®l D + ^W' a ®R^\r) 

a=l a=l 

g g 
+ (2Tr) 2 5(p)5(p') (l\ [V ® l D + ^ W' a ® 1Q \r) (/| [l D ® V + ^ R a ® VQ \r) . 

The momentum eigenstates |$ P (V, {Wq,})) cannot be normalized to unity in the thermody- 
namic limit, but rather satisfy a ^-normalization. For p — p' — 0, there is an additional 
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divergence which is stronger than the ^-normalization. It can be related to the overlap 

between the |$ P (V, {W^})) and the original cMPS {R a })), which is given by 

<? 

{R a }) \%(V, {W a })) = 2n6(p) (l\ \ [V <g> 1 D + ^ W a ® 7Q |r) . (94) 

a=l 

As before, a one-parameter family of local gauge transformations g{x;s) = exp(sh(x)) 
with h(x) G Ql{D; C) induces a map to the kernel of the representation $ p of T$ by setting 
h(x) = he ipx , so that 

\%(M.®\h), {J^{h)}; Q, {R a })) = 0, 

with 

Mf p \h) = [Q, h] + iph and = [Ra, h). (95) 

We henceforth omit the superscript notation of Q and R a . The dimension of the kernel of 
the map $ p is thus D 2 - dimensional, except at p — 0. This can easily be proven, since for 
every non-zero h G gi(D; C), M$ (/i) ^ or 3\f ai $ p (/i) 7^ 0, Va = 1, . . . , N. Indeed, suppose 
that M$ p (/i) = and N<s> P (h) = 0. Imposing that 

N 

0=1 

results in T\hr) = \p\hr) which has no non-trivial solution except at p — 0, where we 
find h = ct-r> with c G C. At nonzero momenta, we can use a gauge fixing condition 
to reduce the number of parameters by D 2 . At p = 0, we can only reduce the number 
of parameters by D 2 — 1 through gauge fixing. But imposing orthogonality to \ty(Q,R)) 
manually at p = allows to discard one additional parameter. For any momentum p, we 
can uniquely fix the gauge of any tangent vector in by setting (l\V ® Id + W ® R = 
or V ® Id + W ® R\r) = 0, corresponding to the left and right gauge fixing conditions 
respectively. It can indeed be checked that with either one of these conditions being satisfied, 
the overlap (&(Q, {R a })\$ p (V, {W a })) given in Eq. ( )94|) vanishes even for p = 0. In addition, 
if either gauge fixing condition is satisfied, the overlap between two tangent vectors simplifies 
significantly, as only the local term survives. Also note the difference with the approach for 
translation non-invariant systems in the previous subsection. There we could impose the 
left or right gauge fixing condition for any x, without this automatically implying that 
\$[V, {W a } , W-&]) _L \^[Q, {R a }}), since a non-zero overlap between the tangent vector and 
the original cMPS could be encoded in the changing boundary vector iu R . 
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IX. CONCLUSION AND OUTLOOK 



This manuscript provides a detailed description of a variational class of wave functions 
for one-dimensional quantum field theories, that goes by the name of "continuous matrix 
product states" . We reviewed different alternative constructions that produce the same class 
of states and have their own merits, e.g. in offering clear hints on how to generalize this 
class to different settings such as open quantum systems or higher-dimensional theories. 

We illustrated how to formulate the cMPS ansatz for the most general class of theories 
including an arbitrary number of bosonic and fermionic particles, and were naturally led to 
a set of constraints that the variational parameters needed to satisfy in order to produce a 
finite kinetic energy density. We also discussed other physical constraints such as fermion 
parity. We then proceeded by explaining in detail how to compute expectation values, 
in particular for the case of systems with open boundary conditions. We provided some 
additional details for the case of systems with translation invariance, where we can use the 
expectation value of a correlation function to define an ultraviolet cutoff within the cMPS 
state. 

We also discussed the important topic of gauge invariance in the cMPS representation. 
Finally we introduced the concept of cMPS tangent vectors, and discussed how the gauge 
invariance allows to represent them in such a way that the metric of the cMPS manifold 
simplifies tremendously. 

While we have not introduced any practical algorithms or recipes for finding cMPS ap- 
proximations of ground states or for describing other physical phenomena, we have intro- 
duced all necessary definitions and concepts in order to comfortably work with cMPS. This 
set of definitions can now be used in follow-up papers that will focus on new algorithms. As 
such, the current paper provides a stepping stone that will hopefully spur more research in 
the context of variational methods for quantum field theories in one dimension and beyond. 
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Appendix A: A useful formula 



Consider an operator U (x, y) denned as 

U(x, y) = y exp 



A(z) dz 



where A is not necessarily antihermitian. This operator satisfies 



— U(x,y) = -A(x)U(x,y), 



—U(x,y) = +U{x,y)A(y). 
dy 



For the derivatives of the inverse operator U(x,y) 1 we can use the general result 

^-U^y)- 1 = -U&y)- 1 (^-U(x,y)) U{x,y)~ l = +U(x, y)^A{x), 
dx \ dx J 

^U&V)- 1 = -Uix/y)- 1 {^U(x,y^j U&y)- 1 = -A(y)U(x,y)-\ 



(Al) 

(A2) 

(A3) 
(A4) 
(A5) 



Now define the following operator quantity depending on an arbitrary operator B 

C(x,y) = U(x,y)BU(x, y y 1 . (A6) 

By taking the derivative with respect to y, we obtain 

^C(x,y) = U(x,y) \A{y),B] U(x,y)-\ 
dy L J 

Integrating dC(x,z)/dz for z from x to y and making use of the initial value C(x,x) = B 

results in 



C(x,y) = B+ / U(x,z) A(z),B U(x,z)~ 1 dz. 



(A7) 

We then multiply this equality with U (x, y) to the right and make use of the obvious identity 
U (x, y) = U(x, z)U(z, y) for any x < z < y in the integral of the right hand side in order to 
obtain our final result 



U(x,y),B 



U(x,z) A(z),B U(z,y)dz. 



(A* 
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We can further generalize this result. Suppose we have two operators U±(x, y) defined as 



U±(x,y) = Texp 



A 1 (z)±A 2 (z)\ dz 



for arbitrary A\^{z). If we consider the quantity 



(A9) 



(A10) 



then we obtain 



— C(x,y) = U-(x,y) - [A 2 (y),Bjj U(x,y)+\ 

using a similar derivation. Continuing along the same line results in 



BU+{x,y)-U^(x,y)B = / U_(x,x 



B,AAz) 



+ {B,A 2 (z)\)U + (z,y)dz. (All) 



Appendix B: Higher order regularity conditions 

In this appendix we derive additional regularity conditions by considering higher deriva- 
tives of the field operators acting on the ground state. Throughout this appendix, we assume 
that Eq. (125]) is fulfilled and R a (x) has well-behaved higher order derivatives. We now con- 
sider the state (d 2 i/j a (x)/dx 2 ) \^f[Q, {Rp}]), which contains a contribution with infinite norm 
unless 

~dR a 



dx 



(x) + [Q(x),R a (x)],R p (x) 



0. 



(Bl) 



J =F 



where [•, -} T is a commutator (— ) or anticommutator (+) for r) a ^ = ±1. If Q and R a 
obey all equations to have a 'well defined' derivative up to order n, so that the state 
(d n ip(x)/dx n ) \ty[Q, {Rp}}) is normalizable, the sufficient condition to eliminate all harm- 
ful contributions from (d n+1 ip(x) /dx n+1 ) \^f[Q,{Rp}]} is 



in— 1 



m-2 



d^ Ra{x) + d^ [g(x) ' Ra{x)] + dx^ [Q{x) > ^'^^ 

+ ... + [Q(x),[...,[Q(x),R(x)]]..],Rp(x) 



0. (B2) 



We can also impose regularity of the mixed derivatives of the A-particle wave function, 
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by first evaluating x(: a (x)ipp(y) \^[Q, {-R 7 }]) 



Mx)My)mQ,{R,}]) = 

e(y-x)tr \BU a fi{—L/ 2 , x)r}p ia R a {x) Up (x, y) Rp (y) U (y, + V 2 )] l^> 

+ 6(x - y) tr \BU aifj (-L/2, y)Rp(y)U a (y, x)R a (x)U(x, +L/2) 



where a new set of operators U a ,p(x, y) (a, ft — 1, . . . , q) was introduced as 



U a ,p(x,y) = Texp 



dz ^ Q(z) ® i + r] a ^r]p >y R y (z) <g> 
7=1 



(B3) 



Note that the regularity condition in Eq. (125]) is sufficient for the annihilation of two particles 
i/j a (x)i/jp(y) \ty[Q, {i? 7 }]) to be continuous at x — y. By first differentiating with respect to 
x, we obtain 



dj) a 
dx 



x) )Mv)MQ,{Rr}]) 



9(y — x) tr 



BU a ,p(-L/2,x)r]p, t 



dR a 
dx 



[x 



+ [Q{x),R c 



, x 



xUp(x,y)Rp(y)U(y,+L/2) 



|0> 



+ 9(x — y) tr 



BU a ,p(-L/2,y)Rp(y)U a (y,x) 



x (^r(z) + ^(^]) U(x, +L/2) 



|n), 



where we have assumed the regularity condition in Eq. (1251) to hold. This allows one to 
eliminate the fixed insertion of particles at position x as well as the terms obtained from 
differentiating the Heaviside functions {i.e. the terms proportional to 5(x — y)). Such terms 
would indeed arise if ipa(x)il)p{i)) \^f[Q, were not continuous at x = y. If we now also 

differentiate with respect to y, we obtain a divergent contribution 



-5(x — y) tr 



BW a ,f,(-L/2,x) 



Rp(x),^(x) + [Q(x),R a (x)} 



U(x, +L/2) 



|n>. 



If we differentiated with respect to y first, and then to x, the divergent contribution is 



5(x — y) tr 



BW a ,p{—L/2, x) 



dRp 
dx 



x) + [Q(x),R p (x)],R a (x) 



U{x, +L/2) 



|n) 
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Since we are working under assumption of the regularity condition [Rp(x), R a (x)] T = 
[Eq. (125])]. it is easy to show that [Rp(x),dR a (x)/dx] T = — [dRp(x)/dx, R a (x)]^ and also 
[Rp(x), [Q(x), R a (x)]] T = —[[Q(x), Rp(x)], R a (x)]^, so that both diverging contributions are 
equal. By imposing 

~dR„ 



dx 



x 



[Q(x),R p (x)],R a { 



x) 



T 



d Ft 

Rp(x),^(x) + [Q(x),R a (x)} 



(B4) 



J =F 



the mixed derivative (dijj a (x) /dx)(dijjp(y) /dy) \^[Q{x), {-R 7 }]) is well defined and normaliz- 
able. Note that Eq. ( 1B4|) is identical to Eq. ( 1B1|) . so that regularity of the mixed product 
of two first order derivatives is guaranteed if the second order derivative is regular, or vice 
versa. 

The higher order regularity conditions derived in this appendix put very strong constraints 
on Q and R a that might be hard to satisfy with finite-dimensional matrices. As mentioned 
in the main text, satisfying the original condition in Eq. (125]) . as imposed by the finiteness 
of the kinetic energy, should be sufficient for most practical applications. 
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